########################################
###Liberalization/Regime Change Graph###
########################################

#Figure A.3

#############	
#############	
### t – 1 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.005,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 1)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction+lag.lncapdist+lag.RCcost_sanction:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc1.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 2 ###
#############
#############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.005,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 2)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction2", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction2+lag.lncapdist+lag.RCcost_sanction2:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc2.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction2", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 3 ###
#############
#############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.005,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 3)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction3", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction3+lag.lncapdist+lag.RCcost_sanction3:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc3.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction3", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 4 ###
#############
#############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.05, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.005,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 4)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction4", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction4+lag.lncapdist+lag.RCcost_sanction4:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc4.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction4", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 5 ###
#############
#############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.005,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 5)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction5", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction5+lag.lncapdist+lag.RCcost_sanction5:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc5.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction5", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 6 ###
#############
#############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.005,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 6)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction6", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction6+lag.lncapdist+lag.RCcost_sanction6:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc6.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction6", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 7 ###
#############
#############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.005,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 7)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction7", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction7+lag.lncapdist+lag.RCcost_sanction7:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc7.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction7", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 8 ###
#############
#############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 8)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction8", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction8+lag.lncapdist+lag.RCcost_sanction8:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc8.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction8", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 9 ###
#############
#############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.005,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 9)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction9", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction9+lag.lncapdist+lag.RCcost_sanction9:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc9.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction9", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


##############
##############
### t – 10 ###
##############
##############
########################################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.007,.01)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 10)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.RCcost_sanction10", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.RCcost_sanction10+lag.lncapdist+lag.RCcost_sanction10:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("RCsanc10.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.RCcost_sanction10", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

